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Abstract. We have coded a Boltzmann solver based on 
a finite difference scheme (Sat method) aiming at calcu- 
lations of neutrino transport in type II supernovae. Close 
comparison between the Boltzmann solver and a Monte 
Carlo transport code has been made for realistic atmo- 
spheres of post bounce core models under the assumption 
of a static background. We have also investigated in de- 
tail the dependence of the results on the numbers of radial, 
angular, and energy grid points and the way to discretize 
the spatial advection term which is used in the Boltzmann 
solver. A general relativistic calculation has been done 
for one of the models. We find overall good agreement 
between the two methods. This gives credibility to both 
methods which are based on completely different formu- 
lations. In particular, the number and energy fluxes and 
the mean energies of the neutrinos show remarkably good 
agreement, because these quantities are determined in a 
region where the angular distribution of the neutrinos is 
nearly isotropic and they are essentially frozen in later on. 
On the other hand, because of a relatively small number of 
angular grid points (which is inevitable due to limitations 
of the computation time) the Boltzmann solver tends to 
underestimate the flux factor and the Eddington factor 
outside the (mean) "neutrinosphere" where the angular 
distribution of the neutrinos becomes highly anisotropic. 
As a result, the neutrino number density is overestimated 
in this region. This fact suggests that one has to be cau- 
tious in applying the Boltzmann solver to a calculation 
of the neutrino heating in the hot-bubble region because 
it might tend to overestimate the local energy deposition 
rate. A comparison shows that this trend is opposite to 
the results obtained with a multi-group flux-limited diffu- 
sion approximation of neutrino transport, employing three 
different flux limiters, all of which lead to an underesti- 
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mation of the hot-bubble heating. The accuracy of the 
Boltzmann solver can be considerably improved by using 
a variable angular mesh to increase the angular resolution 
in the semi-transparent regime. 

Key words: radiative transfer - methods: numerical - 
stars: neutron - supernovae: general - elementary parti- 
cles: neutrinos 



1. Introduction 

Neutrino energy transfer to the matter adjacent to the 
nascent neutron star is supposed to trigger the explosion 
of a massive star (M > 8M Q ) as a type II supernova. 
Since the energy released in neutrinos by the collapsed 
stellar iron core is more than 100 times larger than the ki- 
netic energy of the explosion, only a small fraction of the 
neutrino energy is sufficient to expel the mantle and en- 
velope of the progenitor star. Numerical simulations have 
demonstrated the viability of this neutrino-driven explo- 
sion mechanism (Wilson 1982; Bethe & Wilson 1985; Wil- 
son et al. 1986) but the explosions turned out to be sensi- 
tive to the size of the neutrino luminosities and the neu- 
trino spectra (Janka 1993: Janka & Miiller 1993; Burrows 
& Goshy 1993) both of which determine the power of the 
neutrino energy transfer to the matter outside the average 
neutrinosphere. The rate of energy deposition per nucleon 
via the dominant processes of electron neutrino absorp- 
tion on neutrons and electron antineutrino absorption on 
protons is given by: 
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Here Y n — n n /rib and Y p = n p /rib are the number frac- 
tions of free neutrons and protons, respectively; the nor- 
malization with the baryon density rib indicates that the 
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rate per baryon is calculated in Eq. (1). £^52 denotes the 
luminosity of either v e or v e in units of 10 52 erg/s and is 
the radial position in 10 7 cm. The average of the squared 
neutrino energy, (e 2 15 ), is measured in units of (15 MeV) 2 
and enters through the energy dependence of the neutrino 
and antineutrino absorption cross sections, (fi) is the an- 
gular dilution factor of the neutrino radiation field (the 
"flux factor" , which is equal to the mean value of the co- 
sine of the angle of neutrino propagation relative to the 
radial direction) which varies between values much less 
than unity deep inside the protoneutron star atmosphere, 
about 0.25 around the neutrinosphere, and 1 for radially 
streaming neutrinos very far out. The factor (fi) deter- 
mines the local neutrino energy density according to the 
relation E v = L v /(4nr 2 c (fi)) and thus enters the heating 
rate of Eq. (1). Only far away from the neutrino emitting 
star, (fi) — > 1, and E u dilutes like r~ 2 . 

Although it was found in two-dimensional simulations 
that convective instabilities in the neutrino-heating re- 
gion can help the explosion (Herant et al. 1994; Janka & 
Muller 1993, 1996; Burrows et al. 1995; Miller et al. 1993; 
Shimizu et al. 1994) by the exchange of hot gas from the 
heating layer with cold gas from the postshock region, 
the strength of this convective overturn and its impor- 
tance for the explosion is still a matter of debate (Janka 
& Muller 1996, Mezzacappa et al. 1998, Lichtenstadt et 
al. 1998). In addition, it turns out that the development 
of an explosion remains sensitive to the neutrino lumi- 
nosities and the mean spectral energies even if convec- 
tive overturn lowers the required threshold values. This is 
the case because convective instabilities can develop suffi- 
ciently quickly only when the heating is fast and an unsta- 
ble stratification builds up more quickly than the heated 
matter is advected from the postshock region through 
the gain radius (which is the radius separating neutrino 
cooling inside from neutrino heating outside) down onto 
the neutron star surface (Janka & Muller 1996, Janka & 
Keil 1998). 

"Robust" neutrino-driven explosions might there- 
fore require larger accretion luminosities (to be pre- 
cise: a larger value of the product £^(e 2 ) in Eq. (1)) 
during the early post-bounce phase, or might call 
for enhanced neutrino emission from the core. The 
latter could be caused, for example, by convec- 
tive neutrino transport within the nascent neutron 
star (Burrows 1987; Mayle & Wilson 1988; Wilson & 
Mayle 1988, 1993; Keil et al. 1996) or, alternatively, by 
a suppression of the neutrino opacities at nuclear 
densities through nuclcon correlations (Sawyer 1989; 
Horowitz & Wehrberger 1991; Raffelt & Seckel 1995; Raf- 
felt et al. 1996; Keil et al. 1995; Janka et al. 1996; Bur- 
rows & Sawyer 1998a,b; Reddy et al. 1998a), nuclcon re- 
coil and blocking (Schinder 1990) and/or nuclear in- 
teraction effects in the neutrino-nucleon interactions 
(Prakash et al. 1997; Reddy et al. 1997, 1998b), all of 
which have to date not been taken into account fully self- 



consistently in supernova simulations. The diffusive prop- 
agation of neutrinos out from the very opaque inner core 
is determined by the value of the diffusion constant and 
thus sensitive to these effects. 

Most of the current numerical treatments of neutrino 
transport, however, are deficient not only concerning their 
description of the extremely complex neutrino interac- 
tions in the dense nuclear plasma but also concerning 
their handling of the transition from diffusion to free 
streaming. While the core flux is fixed in the diffusive 
regime, the accretion luminosity as well as the spectra 
of the emitted neutrinos depend on the transport in the 
semitransparent layers around the sphere of last scatter- 
ing. Since neutrino-matter interactions are strongly de- 
pendent on the neutrino energy, neutrinos with different 
energies interact with largely different rates and decou- 
ple in layers with different densities and temperatures. 
The spectral shape of the emergent neutrino flux is there- 
fore different from the thermal spectrum at any partic- 
ular point in the atmosphere. Even more, through the 
factor (fi) in the denominator of Eq. (1) the energy de- 
position rate depends on the angular distribution of the 
neutrinos in the heating region. A quantitatively reli- 
able description of these aspects requires the use of so- 
phisticated transport algorithms which solve the Boltz- 
mann equation instead of approximate methods like flux- 
limited diffusion techniques (Janka 1991a, 1992; Mezza- 
cappa & Bruenn 1993a,b,c; Messer et al. 1998). The de- 
tection of electron antineutrinos from SN 1987A in the 
Kamiokandc II (Hirata et al. 1987) and 1MB laboratories 
(Bionta et al. 1987) and the construction of new, even 
larger neutrino experiments for future supernova neutrino 
measurements have raised additional interest in accurate 
predictions of the detectable neutrino signals from type II 
supernovae. 

Neutrino transport in core collapse supernovae is a 
very complex problem and difficult to treat accurately 
even in the spherically symmetric case. Some of the ma- 
jor difficulties arise from the strong energy dependence 
of the neutrino interactions, the non-conservative and 
anisotropic nature of the scattering processes such as 
neutrino-electron scattering, the non-linearity of the re- 
action kernels through neutrino Fermi blocking, and the 
need to couple neutrino and antineutrino transport for 
the neutrino-pair reactions. Therefore various simplifi- 
cations and approximations have been employed in nu- 
merical simulations of supernova explosions and neutron 
star formation. The so far most widely used approxima- 
tion with a high degree of sophistication is the (multi- 
energy-group) flux-limited diffusion (MGFLD) (Bowers & 
Wilson 1982, Bruenn 1985, Myra et al. 1987, Suzuki 1990, 
Lichtenstadt et al. 1998) where a flux-limiting parame- 
ter is employed in the formulation of the neutrino flux 
to ensure a smooth interpolation between the diffusion 
regime (where the neutrinos are essentially isotropic) and 
the free streaming regime (where the neutrinos move ra- 
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dially outward). Although the limits are accurately repro- 
duced, there is no guarantee that the intermediate regime 
is properly treated. Since in a quasi-stationary situation 
(e.g., for the cooling protoneutron star) the flux and the 
mean energy of the emitted neutrinos are determined in 
the diffusion regime, little change of these is found when 
the flux-limiter is varied (Suzuki 1990) or the transport 
equation is directly solved, e.g., by Monte Carlo calcu- 
lations (Janka 1991a). This, however, is not true when 
the spectral form is considered, because the spectra are 
shaped in the semitransparent surface-near layers. More- 
over, significant differences are also expected for prob- 
lems where the local angular distribution is important in 
the region between the diffusion and free streaming lim- 
its. Due to the factor (fi) appearing in Eq. (1) the hot- 
bubble heating is such a problem, neutrino-antineutrino 
pair annihilation is another problem of this kind. In fact, 
Monte Carlo simulations (Janka & Hillebrandt 1989a,b; 
Janka 1991a; Janka et al. 1992) have shown that all flux- 
limitcrs overestimate the anisotropy of the radiation field 
outside the neutrinosphere, i.e., (ii) = 1 is enforced too 
rapidly (see also Cernohorsky & Bludman 1994). This 
leads to an underestimation of the neutrino heating in the 
hot-bubble region between neutrinosphere and supernova 
shock (Eq. (1)), and sensitivity of the supernova dynam- 
ics to the employed flux-limiting scheme must be expected 
(Messer et al. 1998, Lichtenstadt et al. 1998). 

Modifications of flux-limited diffusion have been 
suggested (Janka 1991a, 1992; Dgani & Janka 1992, 
Cernohorsky & Bludman 1994) by which consider- 
able improvement can be achieved for spherically 
symmetric, static and time-independent backgrounds 
(Smit et al. 1997), but satisfactory performance for the 
general time-dependent and non-stationary case has not 
been demonstrated yet. Therefore the interest turns to- 
wards direct solutions of the Boltzmann equation for neu- 
trino transport, also because the need to check the appli- 
cability of any approximation with more elaborate meth- 
ods remains. Moreover, the rapid increase of the computer 
power and the wish to become independent of ad hoc con- 
straints on generality or accuracy yield a motivation for 
the efforts of several groups (in particular Mezzacappa 
& Bruenn 1993a,b,c and Messer et al. 1998; more recently 
also Burrows 1997) to employ such Boltzmann solvers in 
neutrino-hydrodynamics calculations of supernova explo- 
sions. 

There arc different possibilities to solve the Boltzmann 
equation numerically, one of which is by straightforward 
discretization of spatial, angular, energy, and time vari- 
ables and conversion of the differential equation into a 
finite difference equation which can then be solved for the 
values of the neutrino phase space distribution function at 
the discrete mesh points. Dependent on the number N of 
angular mesh points, this procedure is called Sat method. 
Since solving the equation is computationally very expen- 
sive, there are limitations to the resolution in angle and 



energy space. Therefore tests need to be done whether a 
chosen (and affordable) number of energy and angle grid 
points is sufficient to describe the spectra well and, in 
particular, to reproduce the highly anisotropic neutrino 
distribution outside the neutrinosphere. 

Another, completely different approach to solve the 
Boltzmann equation is the Monte Carlo (MC) method by 
which the probabilistic history of a large number of sam- 
ple neutrinos is followed to simulate the neutrino trans- 
port statistically (Tubbs 1978; Janka 1987, 1991a; Janka 
& Hillebrandt 1989a). In principle, the accuracy of the 
results is only limited by the statistical fluctuations asso- 
ciated with the finite number of sample particles. Since 
the MC transport essentially does not require the use of 
angle and energy grids, it allows one to cope with highly 
anisotropic angular distributions and to treat with high 
accuracy neutrino reactions with an arbitrary degree of 
energy exchange between neutrinos and matter. However, 
the MC method is also computationally very time con- 
suming, in particular if high accuracy on a fine spatial 
grid or at high optical depths is needed. Therefore it is 
not the transport scheme of one's choice for coupling it 
with a hydrodynamics code. 

In the present work, we make use of the advantages of 
the MC method in order to test the accuracy and relia- 
bility of a newly developed neutrino transport code that 
follows the lines of the Sn scheme described by Mezza- 
cappa & Bruenn (1993a,b,c). In particular, we shall test 
the influence of the number of radial, energy, and angular 
mesh points on predicting the spectra and the radial evo- 
lution of the neutrino flux in "realistic" protoneutron star 
atmospheres as found in hydrodynamic simulations of su- 
pernova explosions (Wilson 1988). Since our investigations 
are restricted to static and time- independent backgrounds, 
we concentrate on generic properties of the transport de- 
scription which should also hold for more general situa- 
tions. The radial evolution of the angular distribution of 
the neutrinos is such a property, because it is primarily 
dependent on the profile of the opacity and the geometry 
of the neutrino-decoupling region, but is not very sensitive 
to the details of the temperature and composition in the 
neutron star atmosphere. Finally, good overall agreement 
of the MC and S n results would strengthen the credibility 
of the MC transport with its limited ability to yield high 
spatial resolution. 

The paper is organized as follows. The details of the 
Boltzmann solver and essential information for the MC 
method are given in Sect. 2. In Sect. 3 we describe the 
background models. Section 4 presents the results of our 
comparative calculations, i.e., neutrino spectra, luminosi- 
ties, and Eddington factors. Some of our calculations are 
also compared against results obtained with a MGFLD 
code developed by Suzuki (1994). The dependence of the 
results from the Sn scheme on the energy, angular, and 
radial grid resolution is discussed, too. It is shown that an 
angular mesh that varies with the position in the star can 
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improve the angular resolution and the representation of 
the beamed neutrino distributions without increasing the 
number of angular mesh points. Finally, a summary of our 
results and a discussion of their implications can be found 
in Sect. 5. 

2. Numerical methods 

2.1. Boltzmann solver 
2.1.1. Basic equations 

Our Sat code is based on a finite difference form of the 
general relativistic Boltzmann equation for neutrinos. We 
assume spherical symmetry of the star throughout this pa- 
per. For the Misner-Sharp metric (Misner & Sharp 1964): 
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the Boltzmann equation in the Lagrangian frame takes 
the following form: 
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where p is the cosine of the angle of the neutrino momen- 
tum with respect to the outgoing radial direction and e v 
is the neutrino energy, pb is the baryonic mass density. 

The right hand side of Eq. (3) is the so-called collision 
term, which actually includes absorption, emission, scat- 
tering and pair creation and annihilation of neutrinos, de- 
tails of which are described below. The left hand side, on 
the other hand, looks a little bit different from the form 
used, for example, in Mezzacappa et al. (1993a,b,c). This 
is not only because it is fully general relativistic but also 
because all the velocity dependent terms are expressed as 
time derivatives so that it can easily be coupled to the im- 
plicit general relativistic Lagrangian hydrodynamics code 
(Yamada 1997). In this way a fully coupled, implicit sys- 
tem of the radiation-hydrodynamics equations is formed, 
in which the time derivatives can be treated easily because 
they are just off-diagonal components of the matrix set up 
from the linearized equations. It should be noted, however, 
that since we assume that the matter background is static 
in this paper, all these time derivatives are automatically 
set to be zero, although these terms have been already 
implemented in the code. 

The conserved neutrino number N v in the absence of 
source terms is represented in terms of the chosen inde- 
pendent variables as 



N v = J f v (t, m, n, e v ] 



dm 2i:e 2 de v dp 
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This suggests to cast Eq. (3) in a conservation form 
with respect to the neutrino number. It is also evident 

that the combination of [ — J is more convenient to be 

\PbJ 

used than /„ itself. In the following, therefore, we define 
the specific neutrino distribution function as in Mezza- 
cappa et al. (1993a,b,c) by 

'«-(£) (5) 

and use this quantity as the dependent variable to be 
solved for. 



In the above formula c and G are the velocity of light and 
the gravitational constant, respectively, t is the coordi- 
nate time and m is the baryonic mass coordinate which is 
related to the circumference radius r through the conser- 
vation law of the baryonic mass. In view of combination 
with a Lagrangian hydrodynamics code (Yamada 1997), 
the baryonic mass is chosen to be the independent vari- 
able instead of the radius. 4>, A and r are the metric com- 
ponents which are determined by the Einstein equations. 
In this paper, however, these quantities are given from the 
background models and set to be constant with time dur- 
ing the neutrino transport calculations. f u is the neutrino 
phase space distribution function. Under the assumption 
of spherical symmetry, /„ is a function of t, m, p and e v , 



2.1.2. Finite difference scheme of Boltzmann equations 

As mentioned above the specific neutrino distribution 
function F v is a function of t, m, p and e v . This four- 
dimensional phase space is discretized and Eq. (3) is writ- 
ten as a finite-difference equation. In the time direction we 
adopt a fully implicit differencing. The discretized specific 
distribution function F™- k is defined at the mesh centers 
of the spatial, angular and energy grids. Here the sub- 
scripts k refer to the spatial, angular and energy grid 
points, respectively. The superscript n corresponds to the 
time step. The value at each cell interface is evaluated 
by interpolation of the distribution at two adjacent mesh 
centers. 
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Our finite difference method is essentially the same 
as that of Mezzacappa et al. (1993a,b,c) with some mod- 
ifications. For the spatial advection, the upwind differ- 
ence and the centered difference are linearly averaged with 
the weights determined by the ratio of the mean free 
path to the distance to the stellar surface unlike Mez- 
zacappa et al. (1993a,b,c) who used the ratio of the mean 
free path to the local mesh width. In fact, in the latter 
case we found that the upwind distribution was given too 
large a weight in the optically thick region and thus the 
flux was overestimated. This issue will be addressed later. 

The angular mesh is determined such that each mesh 
center and cell width correspond to the abscissas and 
weights of the Gauss-Legendre quadrature, respectively. 
In angular direction, the neutrino distribution at each in- 
terface is simply taken as the upwind value. The advection 
in the energy space is also approximated by an upwind 
scheme following Mezzacappa & Matzner (1989), although 
this does not allow to conserve both lepton number and 
energy in non-static situations (which are not considered 
here) unless a large number of energy zones is used (Mez- 
zacappa, private communication). 

In typical calculations, 105, 6 and 12 mesh points are 
used for the spatial, angular and energy discretizations, 
respectively. The dependence of the results on the num- 
bers of mesh points will be discussed below. The finite- 
diffcrcnccd Boltzmann equation forms a nonlinear coupled 
system of equations for all radial grid points which is lin- 
earized and solved iteratively by using a Newton-Raphson 
scheme. The linearized equations adopt a block tridiag- 
onal matrix form, which can be efficiently solved by the 
Feautrier method. 

2.2. Monte Carlo method 

Different from the finite difference method, the Monte 
Carlo method constructs the statistical ensemble average 
by following the destinies of individual test particles and 
performing the average when all particles have been trans- 
ported. Due to the fact that neutrinos are fermions it 
is impossible to propagate them independently. Instead, 
the full time-dependent problem has to be simulated by 
following a large number (typically <~ 500,000) of sam- 
ple particles along their trajectories simultaneously in or- 
der to be able to construct the local phase space occupa- 
tion functions and to include anisotropies as well as phase 
space blocking effects self-consistently into the calculation 
of the reaction rates and source terms. The modeling of 
the phase space distribution function from the local parti- 
cle sample must guarantee the correct approach to chem- 
ical equilibrium. Also the Pauli exclusion principle has to 
be satisfied by the statistical average. We refer readers to 
Janka (1987) and references therein for details. 

The stellar background is divided into 15 equispaced 
spherical shells of homogeneous composition and uniform 
thermodynamical conditions, the number of which was 



determined both from physical requirements for spatial 
resolution and from the requirement to have acceptably 
small statistical errors in the local neutrino phase space 
distributions constructed from the chosen number of sam- 
ple particles. Although the Monte Carlo method is essen- 
tially mesh free, about 60 energy bins and approximately 
35 angular bins are used only for representing the phase 
space distribution functions and for calculating the reac- 
tion kernels. Neutrinos are injected into the computational 
volume at the inner boundary in the way described in the 
next section, while particles passing inwardly through the 
inner boundary are simply forgotten. The outer bound- 
ary is treated as a free boundary, where particles escape 
unhindered and no neutrino is assumed to come in from 
outside. 

2.3. Boundary conditions 

In the presented calculations we are mainly interested in 
the neutrino transport in the region where the neutri- 
nos decouple from the stellar background and the emitted 
spectra form. Therefore we calculate the neutrino trans- 
port only in the vicinity of the "neutrinosphere" . For this 
reason we have to set an inner boundary condition as well 
as an outer boundary condition for each model. At the 
outer boundary we impose the condition that no neutri- 
nos enter the computational volume from outside. In the 
Boltzmann solver, this is realized by setting 



f , „ e x J/4w, Mi £«/) for fi > 
M r s , it, e v ) - | Q for ^ < ^ 



(6) 



and r s is the radius of the outer surface which is dislocated 
outward from the outermost mesh center by half a radial 
cell width. 

On the other hand, we have to specify the distribu- 
tion of the neutrinos coming into the computational vol- 
ume at the inner boundary ne which is dislocated in- 
ward from the innermost mesh center r m ; n by half a ra- 
dial cell width. For this purpose we adopt Fermi-Dirac 
distribution functions, in which the temperature, chemi- 
cal potential and a normalization factor are determined 
such that the neutrino number density at hb (where the 
neutrinos are essentially isotropically distributed because 
the inner boundary is chosen to be located at high opti- 
cal depth), the average energy (ej) and the width of the 

(4) 



energy spectrum, measured by the parameter 



(see 



Janka & Hillebrandt 1989b), are reproduced as given by 
Wilson's (1988) models. Thus the inner boundary condi- 
tion is set as : 
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where A, /^ib and Tib are the fitting parameters, the values 
of which are summarized in Table 1 for all considered mod- 
els and neutrino species. Concerning the distribution of 
neutrinos that leave the computational volume, we impose 
in Eq. (7) in the Boltzmann solver the condition that it is 
the same as the distribution at the innermost mesh cen- 
ter. From the physical point of view, however, and treated 
correctly in the Monte Carlo simulations, it should be de- 
termined by the fraction of neutrinos which is emitted or 
backscattered towards the inner boundary. This is in gen- 
eral different from what the phase space distribution at 
the innermost mesh center yields because the mean free 
path of the neutrinos near the inner boundary is shorter 
than the mesh width. As a result the imposed condition 
for \x < in Eq. (7) leads to a minor discrepancy of the 
treatment of the inner boundary condition in the Monte 
Carlo and Boltzmann computations and sometimes causes 
a small oscillation of the neutrino distribution near the in- 
ner boundary in the latter computations. This issue will 
be revisited later. 



Table 1. Fitting parameters for the Fermi-Dirac distributions 
which describe the spectra at the inner boundary. 



models 


neutrinos 


Tib [MeV] 


(Uib [MeV] 


A 


Wl 




9.56932 


7.85190 


1.221490 






9.38259 


7.13610 


3.874140 






10.32539 


30.87195 


17.911456 


W2 




10.52454 


5.75932 


3.672680 


W3 


Ve 


9.87875 


11.27370 


7.108076 



2.4- Input physics 
2.4.1. Neutrino reactions 

The neutrino opacities of dense neutron star matter are 
still one of the major uncertainties of supernova simula- 
tions. Theoretical and numerical complications arise from 
the description and treatment of nucleon thermal mo- 
tion and recoil (Schinder 1990), nuclear force effects and 
nucleon blocking (Prakash et al. 1997; Reddy et al. 1997), 
and nucleon correlations, spatial (Sawyer 1989; Burrows & 
Sawyer 1998a,b; Reddy et al. 1998a) as well as temporal 
(Raffelt & Seckel 1995; Raffelt et al. 1996). Although in 
particular nucleon recoil and auto-correlations might play 
an important role even in the sub-nuclear outer layers of 
the protoneutron star down to densities below 10 13 gcm~ 3 
(Janka et al. 1996; Hannestad & Raffelt 1997) we do not 
concentrate on this problem here but rather employ the 
standard description of the neutrino opacities, accord- 
ing to which neutrinos interact with isolated nucleons 
(see, e.g., Tubbs & Schramm 1975; Bruenn 1985). Also, 



as in most other simulations, bremsstrahlung production 
of neutrino-antineutrino pairs is neglected here, although 
it may be important as pointed out by Suzuki (1993) and 
more recently by Burrows (1997) and Hannestad & Raf- 
felt (1997). Doing so, we intend to enable comparison with 
other (already published) work and want to avoid the mix- 
ing of effects from a different numerical treatment of the 
transport with those from a non-standard description of 
neutrino- nucleon interactions or from the inclusion of pro- 
cesses typically not considered in the past. 

The following neutrino reactions have been imple- 
mented in our codes. 

[1] v e + n ^ e + p 

(electron- type neutrino absorption on neutrons), 

[2] v e + p # e+ + n 

(electron- type anti- neutrino absorption on protons), 

[3] v + N ^ v + N 

(neutrino scattering on nucleons), 

[4] v + e # v + e 

(neutrino scattering on electrons). 

In addition to the above reactions, the following processes 
have also been implemented in both codes, although these 
reactions are not used in the present paper. 

[5] v e + A # A + e~ 

(electron-type neutrino absorption on nuclei), 

[6] v + A # v + A 

(neutrino coherent scattering on nuclei), 

[7] e" + e+ # v + v 

(electron-positron pair annihilation and creation), 

[8] 7* # v + v 

(plasmon decay and creation). 

Since the pair processes are not taken into account in this 
paper, we can treat each species of neutrinos separately. A 
test showed that in the considered protoneutron star at- 
mospheres neutrino pair creation and annihilation as well 
as processes involving nuclei are unimportant to determine 
the fluxes and spectra. 

2.4.2. EOS 

In this paper we use a simplified equation of state, in which 
only nucleons, electrons, alpha particles and photons are 
included. They are all treated as ideal gases. For given 
density, temperature and electron fraction we derive the 
mass fractions and the chemical potentials of nucleons and 
the electron chemical potential from this equation of state. 
The disregard of nuclei is well justified for the densities 
and temperatures we are considering, where most of the 
nuclei are dissociated into free nucleons. This is actually 
confirmed by comparing our EOS with the more realis- 
tic EOS of Wolff that is based on the Skyrme-Hartree- 
Fock method (Hillcbrandt & Wolff 1985). Only very small 
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differences of the nucleon chemical potentials are found 
for the innermost region where small amounts of nuclei 
appear and for the outermost region where some con- 
tribution from alpha particles is mixed into the stellar 
medium. We also repeated some of the calculations mak- 
ing use of the Wolff EOS with the nuclei-related reactions 
[5], [6] (neutrino emission, absorption and scattering on 
nuclei) turned on and found qualitatively and quantita- 
tively the same results. Hence the nuclei-related reactions 
are switched off in the calculations described below. 

3. Models 

3.1. Stellar models 

The time-dependent transport calculations presented here 
were performed for background profiles which are rep- 
resentative of protoneutron star atmospheres during the 
quasi-static neutrino cooling phase (Wilson 1988). At this 
stage, several seconds after core bounce, the typical evolu- 
tion timescale of density, temperature, and electron frac- 
tion is much longer than the timescale for neutrinos to 
reach a stationary state. Therefore our assumption of a 
static and time-independent background is justified. In 
addition, our interest is focused on the radial evolution 
of the Eddington factors and on a test of the influence of 
the energy and angle resolution used in the S n Boltzmann 
solver. Both aims do not require a fully self-consistent ap- 
proach which takes into account the evolution of the stellar 
background (in particular of the temperature and compo- 
sition). In fact, the Eddington factors are normalized an- 
gular integrals of the radiation intensity and as such reflect 
very general characteristics of the geometrical structure of 
the atmosphere where neutrinos and matter decouple. 

Profiles from Wilson's (1988) protoneutron star model 
were taken for three different times, 3.32, 5.77, and 7.81 s 
after core bounce. With the chosen fundamental variables 
density p, temperature T, and electron fraction Y e , the 
thermodynamical state is defined for the plasma consist- 
ing of non-relativistic free nucleons, arbitrarily relativistic 
and degenerate electrons and positrons, and photons in 
thermal equilibrium. Figures 1-3 show the input used for 
the three models. In Fig. 3 also the general relativistic 
metric coefficients y/gu — and ^—g r r — e A are given 
as provided by Wilson's data and used for a comparative 
general relativistic calculation of the neutrino transport 
in model W3. 

3.2. Computed transport models 

All models computed with the Boltzmann code are sum- 
marized in Table 2. Models ST are the standard models, in 
which 105 uniform spatial, 6 angular and 12 energy mesh 
points were used. The energy mesh is logarithmically uni- 
form and covers 0.9 - 110 MeV. The numbers of angu- 
lar grid points and energy grid points were increased in 
models FA and FE, respectively. In model CS we used the 




1.25 1.3 1.35 1.4 1.45 1.5 1.55 1.6 
r Ll0 6 cmj 



Fig. 1. The profile of density, temperature and electron 
fraction for Wilson's (1988) post bounce core model Wl 
(t = 3.32 s). 
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Fig. 2. The profile of density, temperature and electron frac- 
tion for Wilson's post bounce core model W2 (t = 5.77 s). 
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Fig. 3. The profile of density, temperature and electron frac- 
tion as well as the metric coefficients yfgtt (open squares) and 
a/— g rr (filled squares) for Wilson's post bounce core model W3 
(t = 7.81s). 



same radial grid as in the Monte Carlo simulations where 
15 radial zones were chosen. 105 spatial mesh points were 
again used in model NI with no interpolations of density, 
temperature and electron fraction in the radial grid of the 
Monte Carlo simulations. Model GR took into account the 
general relativistic effects. We used a non-uniform spatial 
mesh in model NU. A different interpolation of up-wind 
differencing and centered differencing was tried for the ra- 
dial advection term in model DI. We assumed the nucleon 
scattering to be isotropic in model IS. As is understood 
from Table 2, most of the comparisons were done for the 
electron-type anti-neutrinos, since they are most impor- 
tant from the observational point of view. 



4. Numerical results 

4-1. Luminosity and average energy 

In the following, we consider the neutrino transport results 
after the time-dependent simulations have reached steady 
states. 

First we compare observable quantities such as the lu- 
minosity, average energy and average squared energy of 
the neutrino flux for models ST and GR. These quantities 



are calculated at the outermost spatial zone as 
27re 2 ,de„d/i 



L u = 



47rr 2 c 



I 



{e v ) = 



(4) = 



(he) 3 



(he) 



(hep 
27T£^de„d/i 



fv M 



(he) 



f, 



li el 



I 



2i:e v de v dii 
(he] 3 



(8) 



(9) 



(10) 



where r s is again the surface radius, and the redshift cor- 
rections from the surface up to infinity are not taken into 
account. The results are summarized in Table 3. 

The electron-type neutrino has the lowest energy while 
the muon-type neutrino has the highest. The reason for 
this is that the electron-type neutrino has the shortest 
mean free path due to absorptions on the abundant neu- 
trons, and decouples in the surface-near layers where the 
temperature is lower. In contrast, the muon and tau neu- 
trinos do not interact with particles of the stellar medium 
by charged currents and therefore their thermal decou- 
pling occurs at a higher temperature. The luminosity for 
the electron-type anti-neutrino gets smaller as time passes, 
which is due to the cooling and shrinking of the protoneu- 
tron star. As can be seen in the table, the agreement of all 
quantities between the two methods is very good, which 
confirms the statistical convergence of the Monte Carlo 
simulations. The average energy is in general determined 
accurately because the energy spectrum is shaped in the 
region where the neutrino angular distribution is not very 
anisotropic. Moreover, possible effects due to the rather 
coarse angular resolution of the Boltzmann code essen- 
tially cancel out by taking the ratios of Eqs. (9) and (10). 
On the other hand, the luminosity and the number flux of 
the Monte Carlo computations are also well reproduced by 
the Boltzmann results. This is due to the fact that these 
quantities are also determined deep inside the star and 
are nearly conserved farther out. In fact, when integrating 
Eq. (3) over angle and energy multiplying with unity and 
e„, respectively, ignoring all time derivatives and general 
relativistic effects, one gets 



d m 
dLjjr) 
d m 



1 

Pb 



Pb J 



27T£ 2 ,d£„d/u 

(he) 5 
27re^de l/ d/u 

(he) 5 



St 



coll 
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St 



(11) 
(12) 



coll 



where F™{r) and L e v (r) are the number and energy fluxes 
of neutrinos at radius r, respectively, and defined as, 

27r£ld£ ^.Up , (13) 
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Table 2. Characteristics of the calculated models. 105 r x 6 M x 12 £ implies that 105 spatial mesh points, 6 angular mesh points 
and 12 energy mesh points are used. See the text for details. 



model 


mesh 


background model 


v species 


notes 


ST1 


105 r x 6 M x 12 £ 


Wl 


v e 


interpolation in spatial grid of MC 


ST2 


105 r x 6 M x 12 e 


Wl 


Ve 


interpolation in spatial grid of MC 


ST3 


105 r x 6 M x 12 e 


W2 


V e 


interpolation in spatial grid of MC 


ST4 


105 r x 6 M x 12 e 


W3 


Ve 


interpolation in spatial grid of MC 


ST5 


105 r x 6 M x 12 e 


Wl 


Vy. 


interpolation in spatial grid of MC 


FA 


105 r x 10 M x 12 e 


Wl 


Ve 


interpolation in spatial grid of MC 


FE 


105 r x 6 M x 18 e 


Wl 


V e 


interpolation in spatial grid of MC 


CS 


15 r x 6 M x 12 e 


Wl 


v e 


same spatial grid as in MC 


NI 


105 r x 6 M x 12 e 


Wl 


V e 


no interpolation in spatial grid of MC 


GR 


105 r x 6 M x 12 e 


W3 


V e 


general relativity included 


NU 


105 r x 6 M x 12 e 


Wl 


V e 


non-uniform spatial mesh 


DI 


105 r x 6 M x 12 e 


Wl 




different interpolation for conservative radial advection 


IS 


105 r x 6 M x 12 e 


Wl 


v e 


isotropic v-N scattering 



Table 3. Luminosity, average energy and average squared energy for models ST1-ST5 and for model GR. B and MC in the 
second column refer to the Boltzmann simulations and the Monte Carlo simulations, respectively. For the definitions of L v , (e v ) 
and (el^ see Eqs. (8), (9) and (10) in the text. 







Ve 


Ve 


V £ 


Ve 


VfX 


Ve 






Wl 


Wl 


W2 


W3 


Wl 


GR, W3 




B 


7.1 x 10 51 


4.0 x 10 51 


2.5 x 10 51 


1.5 x 10 51 


4.9 x 10 51 


1.3 x 10 51 


[erg/sec] 


MC 


7.3 x 10 51 


4.0 x 10 51 


2.6 x 10 51 


1.5 x 10 51 


4.9 x 10 51 


1.3 x 10 51 




B 


12.8 


16.3 


15.9 


16.2 


24.3 


15.5 


[MeV] 


MC 


12.7 


16.1 


15.8 


16.3 


24.2 


15.6 


[MoV 2 ] 


B 


198.5 


329.0 


309.9 


324.5 


727.8 


300.3 


MC 


198.6 


322.6 


308.9 


327.1 


724.3 


300.6 



As expected intuitively, the scattering kernels drop out of 
the right hand side of Eq. (11) while only the isoenergetic 
scattering does not contribute on the right hand side of 
Eq. (12). In the given Boltzmann code Eqs. (11) and (12) 
are discretizcd in a conservative form for the radial ad- 
vection. Therefore it is clear that the Boltzmann code can 
calculate the number and energy fluxes accurately if also 
the number and energy exchange by the reactions are cal- 
culated accurately in the source terms on the right hand 
sides of the equations. Moreover, the radial evolutions of 
F™(r) and L e v (r) are entirely determined by the number 
and energy exchange through reactions of which the net 
effect is small in an atmospheric layer which is in a station- 
ary state and reemits as much energy and lepton number 
as it absorbs. Changes of the number fluxes and lumi- 
nosities in the considered protoneutron star atmospheres 
occur only in regions where the neutrino distribution is 
still essentially isotropic and possible effects due to an in- 
sufficient angular resolution in the Boltzmann Sjv scheme 
do not cause problems. For all these reasons it is not sur- 
prising that the same quality of agreement is also found 



for the radial evolutions of the luminosity, average energy 
and average squared energy and that this is also true for 
the general relativistic case. 

4-2. Neutrino energy spectra 

In the previous section we discussed only the energy 
and angle integrated quantities. However, we also pro- 
vide information about the energy spectra of each neu- 
trino species, because they yield more evidence about the 
quality of the agreement between the calculations with the 
different codes. 

Figs. 4 and 5 show energy flux spectra defined as 

1 dL v _ f 2ndfi 3 

i^&^y (h^ Ufi£u ' ( j 

at the protoneutron star surface for different cases. For 
the reasons discussed in section 4.1, the spectra computed 
with the Boltzmann code (symbols) and the Monte Carlo 
code (lines) show excellent agreement. The number and 
the distribution of the energy bins in the Boltzmann code 
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seem to be adequate to reproduce the highly resolved 
Monte Carlo spectra. 
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Fig. 4. The neutrino energy flux spectra at the protoneutron 
star surface. The lines show the results of the Monte Carlo 
simulations (the thick solid line, the short dashed line and 
the thin solid line for electron-type neutrinos, electron-type 
anti-neutrinos and muon-type neutrinos, respectively, in case 
of background model Wl, and the dotted line and the long 
dashed line for electron-type anti-neutrinos in case of back- 
ground models W2 and W3, respectively). The symbols rep- 
resent the results of the corresponding Boltzmann simulations 
(the filled triangles, the filled circles and the open triangles 
for electron-type neutrinos, electron-type anti-neutrinos and 
muon-type neutrinos, respectively, in case of background model 
Wl, and the filled squares and the filled diamonds for elec- 
tron-type anti-neutrinos in case of background models W2 and 
W3, respectively). 



4-3. Flux factor and Eddington factor 

So far we discussed only angle integrated quantities since 
they are observable. However we are also interested in the 
angular distributions of the neutrinos, because informa- 
tion about the angular distributions is important to deter- 
mine the neutrino heating rate in the hot-bubble region 
(see Eq. (1)). Although it is an advantage of the Boltz- 
mann solver over MGFLD that one does not have to as- 
sume an ad hoc closure relation between the angular mo- 
ments of the distribution function, one should remember 
that the usable number of angular mesh points is severely 



0.5 



- 
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0.2 0.4 0.6 0.8 1 

neutrino energy | 1 B MeV] 



Fig. 5. The electron-type anti-neutrino energy flux spectra 
for background model W3 with and without general relativ- 
ity. The long dashed line and the filled squares show the re- 
sults of the Monte Carlo simulation and the Boltzmann simula- 
tion, respectively, with the general relativistic effects included. 
The solid line and the filled triangles are for the correspond- 
ing non-relativistic results. The short dashed line is the red- 
shift-corrected energy spectrum at infinity calculated from the 
Monte Carlo results, the other results are given at the pro- 
toneutron star surface. 



limited. In the Fcautricr method the computation time in- 
creases in proportion to the third power of the dimension 
of the blocks in the tridiagonal block matrix which has 
to be inverted when one chooses the radius as the outer- 
most variable of the do- loops. The dimension of one block, 
in turn, is linearly proportional to the number of angular 
mesh points. The same dependence holds for the number 
of energy mesh points and the number of neutrino species. 
In the standard calculations we use 6 angular mesh points 
and 12 energy grid points, and we treat a single neutrino 
species at a time, which corresponds to a block matrix 
size of 72. On the other hand, the number of spatial grid 
points is about 100. The CPU time is a few seconds per 
inversion of the whole matrix on a single vector processor 
of a Fujitsu VPP500. Hence, use of more than 10 angu- 
lar mesh points is almost prohibitive for calculations with 
three neutrino species even with a highly parallelized ma- 
trix inversion routine (Sumiyoshi & Ebisuzaki 1998). It is, 
therefore, important to clarify the sensitivity of the accu- 
racy to the angular resolution. 



Shoichi Yamada et al.: Neutrino transport in type II supernovae: Boltzmann solver vs. Monte Carlo method 



11 



For this reason we consider the flux factor (/x) el and 
the Eddington factor (^ 2 ) cl which are defined as : 



Mel 



fv HEv 



2-KS v de v d[i 



I 



: fv 

27re^de J/ d/i 2 



(hc) s 



(he) 2 



2-KE v de v dii 
(hc]3 



(16) 



(17) 



fv &v 



Here the subscript "el" means that the averages are de- 
fined with the weight of energy. In MGFLD these fac- 
tors are related with each other by a closure condi- 
tion which can be derived from the employed flux-limiter 
(Janka 1991a, 1992). For simplicity we discuss here only 
energy integrated quantities as defined above. The fun- 
damental features are similar for the individual energy 
groups. 

In Figs. 6-8, we show the radial evolutions of the flux 
factors and the Eddington factors for all neutrino species 
in case of background model Wl. The upper panels show 
the flux factors and the lower panels the corresponding 
Eddington factors. The solid lines are the results of the 
Boltzmann simulations (having the finer radial resolution) 
and the filled triangles are those of the Monte Carlo simu- 
lations. As can be seen, near the inner boundary the flux 
factors are almost zero while the Eddington factors are 
1/3, which implies that the neutrino angular distribution 
is nearly isotropic, a consequence of the fact that the neu- 
trinos arc in equilibrium with the surrounding matter. As 
we move outward, both factors begin to deviate from these 
values, reflecting the increase of the mean free path and 
a more rapid diffusion. Farther out, the angular moments 
increase monotonically towards unity, the value in the free 
streaming limit, as the angular distribution gets more and 
more forward peaked with increasing distance from the 
source. As can be seen clearly in Figs. 6-8, the Boltzmann 
solver tends to underestimate both angular moments in 
the outermost region, where the neutrino angular distri- 
bution is most forward peaked. This can be explained by 
the insufficient angular resolution, or, to be more specific, 
by the fact that in case of the employed Gauss-Legendre 
quadrature the maximum angle cosine ^ max of the angular 
grid is significantly less than unity, if not a large number 
of angular grid points is used. This is directly confirmed 
by using a larger number of angular mesh points or, in 
particular, a variable angular mesh (see Sect. 4.4). 

The same trend is also present in the general relativis- 
tic case, model GR, as shown in Fig. 9. It turned out that 
the ray bending effect which tends to isotropize the angu- 
lar distribution of the neutrinos is not very important and 
that differences between the Monte Carlo results and the 
Boltzmann results are significantly larger. 



In Fig. 10 we show the flux factor and the Eddington 
factor for model FA where wc employed 10 angular mesh 
points instead of 6. The long dashed lines are for model FA 
and the short dashed lines depict the result of model ST2, 
the corresponding standard model, for comparison. The 
discrepancy between the Boltzmann simulation and the 
Monte Carlo simulation is reduced with the increase of 
the number of angular mesh points. This supports our 
interpretation that the deviation stems entirely from the 
insufficient angular resolution and/or the unfavorable lo- 
cation of the angular grid points for the Gauss-Legendre 
quadrature used in the Boltzmann simulations. Indeed, 
the degree of improvement is consistent with the fact that 
our finite difference scheme is of first order for the angular 
advection, since we always take the upwind differencing 
as explained above. Even with the finer 10-zone angular 
mesh the deviation of both factors from the exact values 
given by the Monte Carlo result is significant. The maxi- 
mum deviations, 1 — ^ max for the flux factor and 1 — /x^ ax 
for the Eddington factor, are approached as one goes far- 
ther out into the optically thin regime. This is visible in 
Fig. 11 which displays the ratios of the Boltzmann to the 
Monte Carlo results for the flux factors and the Eddington 
factors in case of 6 and 10 angular bins and the variable 
angular mesh. (The relatively large discrepancies of the 
flux factors for smaller radii are explained by slight differ- 
ences of the treatment of the inner boundary condition, 
see Sect. 2.3, and by the fact that the flux factor adopts 
very small values in the optically thick region.) It should 
be mentioned that the flux factors calculated in the recent 
paper by Messer et al. (1998) do not converge to unity but 
saturate at a nearly constant lower level (around 0.9) even 
far outside of the neutrinosphere. This reflects the use of 
AWx = 0.93 for the largest ^-bin of the angular grid in 
the 6-point quadrature of the S§ method. 

It is interesting to see that this tendency of the Boltz- 
mann solver is completely opposite to that of MGFLD. 
Janka (1991a, 1992) pointed out that all flux limiters used 
so far tend to overestimate the flux factor and the Ed- 
dington factor in the optically thin region, which im- 
plies that the neutrino angular distribution approaches 
the free streaming limit much too rapidly (see also Messer 
et al. 1998). In order to confirm this statement, transport 
calculations with MGFLD were done for the same mod- 
els with three different flux limiters, which are Bruenn's 
(BR), Levermore & Pomraning's (LP) and Mayle & 
Wilson's (MW). We refer readers to Janka (1992) and 
Suzuki (1994) and references therein for details on the 
flux limiters. We show in Figs. 12 and 13 the flux fac- 
tors and local number densities of the electron-type neu- 
trino and electron-type anti-neutrino for model Wl, re- 
spectively. It is clear that all flux limiters overestimate the 
forward peaking of the angular distributions of the neutri- 
nos in the optically thin region, a trend that holds for all 
neutrino species and is not dependent on the background 
model. The typical deviation of MGFLD results from the 
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Monte Carlo results is much larger than that between the 
Boltzmann solver and the Monte Carlo method. 
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Fig. 6. The flux factor (upper panel) and the Eddington factor 
(lower panel) of the electron-type neutrino for model Wl. The 
filled triangles and the solid line show the results of the Monte 
Carlo simulation and the Boltzmann simulation, respectively. 



From the lower panels of Figs. 12 and 13 we learn that 
the local neutrino number density, which is given by 



' 27r£^d£„d/x 



(18) 



is overestimated in case of the Boltzmann solver (by about 
10%) and underestimated for MGFLD (by approximately 
30%) in the optically thin region. This is understood from 
the fact that the number flux, which is defined as 



F?(r) = 47rr 2 c 



/ 



27T£ 2 d£„d^ 



is related to the local neutrino number density by 
K(r) 



n v {r) 



4irr 2 c (fj) e0 



(19) 



(20) 



(/x) e0 denotes the average angle cosine for the neutrino 
number flux and is calculated from Eq. (16) with a fac- 
tor e v omitted under the integrals in the numerator and 
denominator. Since the number flux is determined deep 
inside the protoneutron star atmosphere where the neu- 
trinos are still nearly isotropic, and conserved farther out, 




1.4 1.5 

r [10 6 cm] 



Fig. 7. The same as Fig. 6 but for the electron-type 
anti-ncutrino. 
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Fig. 8. The same as Fig. 6 but for the muon-type neutrino. 
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Fig. 9. The flux factor (upper panel) and the Eddington factor 
(lower panel) of the electron-type anti-neutrino for model GR 
where the general relativistic effects are taken into account. 
The solid line represents the results of the Boltzmann sim- 
ulation while the solid triangles correspond to those of the 
Monte Carlo simulation. For comparison the corresponding 
non-relativistic results, model ST4, are shown with the solid 
squares and the dashed lines for the Monte Carlo simulation 
and the Boltzmann simulation, respectively. 

it is not affected by problems with a coarse angular reso- 
lution occurring in the optically thin regime. For this rea- 
son, the number and energy fluxes agree well between the 
Monte Carlo method and the Boltzmann solver irrespec- 
tive of the angular resolution as long as the Boltzmann 
solver is based on conservative finite differencing in the 
radial direction. The good agreement of the fluxes is con- 
firmed by Fig. 14, which depicts the radial behavior of the 
number flux in case of models ST1, ST2 and ST5. 

It is now clear from Eq. (20) that an underestimation 
(overestimation) of the flux factor leads to an overestima- 
tion (underestimation ) of the number density, if the flux 
is the same. It should be noted that even in the outermost 
zone of our computational region in the atmosphere of 
a protoneutron star, the neutrino angular distribution is 
not so strongly forward peaked as in the hot-bubble region 
farther out. Hence it must be expected that the errors by 
an over- or underestimation of the neutrino number den- 
sity might be even larger in the hot-bubble region where 
neutrino heating takes place. 

Since the neutrino heating rate is proportional to the 
local neutrino number density (actually: energy density) — 



Fig. 10. The flux factors (upper panel) and the Eddington fac- 
tor (lower panel) of the electron-type anti-neutrino for back- 
ground model Wl with the long dashed lines for model FA 
where 10 angular mesh points are used instead of 6. The short 
dashed lines show for comparison the results of the correspond- 
ing standard model ST2. The solid lines are the result obtained 
with the variable angular mesh method. The Monte Carlo re- 
sult is shown with the solid triangles. 

this is why the inverse of the flux factor appears in 
Eq. (1) — the application of the Boltzmann solver with a 
relatively small number of angular mesh points may lead 
to an overestimation of the neutrino energy deposition in 
the hot bubble for disadvantageous situations, in particu- 
lar when most of the heating occurs in those regions where 
the deviation of the flux factor from the correct value is 
significant. In contrast, all flux limiters used in MGFLD 
underestimate the heating rate significantly. Therefore, 
even for the Boltzmann solver, improvement of the an- 
gular resolution or a redistribution of the angular grid 
points is desirable in order to a priori avoid inaccurate 
evaluation of the heating rate. As already mentioned, an 
increase of the number of angular mesh points is not feasi- 
ble. Choosing a variable angular mesh which adjusts mesh 
point locations in dependence of time and spatial position 
might be a solution. This issue will be discussed in the 
next subsection. 

4-4- Variable angular mesh in the Boltzmann solver 

Here we attempt to improve the angular resolution of the 
Boltzmann solver by redistributing the angular grid points 
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Fig. 11. Same as Fig. 10 but showing the ratios of the Boltz- 
mann to the Monte Carlo results for the flux factor (up- 
per panel) and the Eddington factor (lower panel). Again, 
the Boltzmann calculations with 6 angular bins (standard 
model ST2, short dashed lines), 10 angular bins (model FA, 
long dashed lines) and the variable angular mesh (solid lines) 
are shown. Between the radial grid points of the Monte Carlo 
simulation (whose locations are indicated by solid triangles), 
the Monte Carlo results are interpolated by cubic splines. 



in dependence of time and position so that their density 
is enhanced in the forward direction in the optically thin 
region where the neutrino angular distribution becomes 
strongly forward peaked and the Boltzmann solver under- 
estimates the flux factor and the Eddington factor. This 
requires adding extra angular advection terms in the nu- 
merical scheme which compensate for the motions of the 
angular mesh points. 

We assume that the position of each interface of the 
angular grid is a function of time, baryonic mass and neu- 
trino energy, i.e., p\ = pi(t, m, e v ). Integrating Eq.(3) over 
angular bins then leads to the following additional angular 
advection fluxes at each angular mesh interface I : 



dt 



dm 



r 2 /„ 



(21) 
(22) 



ia(lnp b r 3 ) 2 



- e^ 47T7^ 



dm 



Fig. 12. The flux factors and number densities of the elec- 
tron-type neutrino as obtained by MGFLD with three different 
flux limiters, Bruenn's (BR) with the short dashed lines, Lev- 
ermore & Pomraning's (LP) with the long dashed lines, and 
Mayle & Wilson's (MW) with the dash-dotted lines. The back- 
ground model is Wl. For comparison, the Monte Carlo result 
and the Boltzmann result (model ST1) are also plotted with 
the triangles and the solid lines, respectively. 



1 d (In r) 
c d t 



d Mi 



(23) 



It is easy to understand that Eqs. (21)-(23) originate from 
the variability of the angular mesh points because of the 
differentials of the m with respect to time, mass and neu- 
trino energy. 

Since the neutrino reaction rates are strongly energy 
dependent and so is the neutrino angular distribution, it 
would be desirable to implement the energy dependent an- 
gular mesh according to Eq. (23). In the current prelim- 
inary attempt, however, we installed only Eqs. (21) and 
(22) for simplicity. Incidentally, since Eq. (23) is propor- 
tional to in static background calculations, it is any- 
way negligible for the models considered here. We note 
that the motion of mesh points is not calculated implic- 
itly, that is, the angular mesh points for the next time 
step are determined from the neutrino angular distribu- 
tion at the current time step and are kept fixed during 
the implicit calculation of the transport for the next step. 

In Fig. 10, we show both the flux factor and the Ed- 
dington factor obtained from the computation with the 
variable angular mesh. The improvement is clear from a 
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Fig. 13. The same as Fig. 12 but for the electron-type 
anti-neutrino. 
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Fig. 14. The number fluxes of all types of neutrinos for back- 
ground model Wl. The filled symbols and the solid lines show 
the results of the Monte Carlo simulations and the Boltzmann 
simulations, respectively. 



comparison with the result of model FA which employed 
10 angular mesh points and is also shown in the figure. It 
should be emphasized that increasing the number of angu- 
lar mesh points from 6 to 10 leads to an increase of CPU 
time by a factor of <~4.5, while the additional operations 
for the variable angular mesh imply negligible computa- 
tional load. 

We repeated all Boltzmann calculations with the vari- 
able angular mesh method and found that the same im- 
provement could be achieved for all cases. Our scheme is 
stable at least for the static background models, although 
the stability for dynamical background models remains 
to be tested. Thus we think this method is promising in 
applying the Boltzmann solver to the study of neutrino 
heating in the hot-bubble region of supernovae, although 
there is room for improvement concerning the prescription 
of the motion of the mesh points and the implementation 
of the energy-dependent angular mesh. 

4-5. Spatial and energy resolution in the Boltzmann solver 
and radial advection 

In this section we discuss how the numerical results change 
in dependence of the number of spatial and energy grid 
points, the boundary condition and the treatment of the 
radial advection in the Boltzmann solver. 

In Fig. 15 we show the energy spectrum of electron- 
type anti-neutrinos for model FE, in which we used 18 
energy mesh points, compared with the spectrum for the 
corresponding standard model ST2 which has 12 energy 
zones. No qualitative or quantitative difference is found 
between the two cases. This is also true for the luminosity 
and the angular distribution. Thus we think that about 
15 energy mesh points are sufficient for the calculation of 
the energy spectrum. These results are in agreement with 
previous findings by Mezzacappa & Bruenn (1993a,b,c). 

In models CS and NI we reduced the spatial resolu- 
tion, because the Monte Carlo simulations were done with 
only 15 radial mesh points which were used to represent 
the stellar background on which the reaction kernels were 
evaluated. Another motivation for testing the sensitivity 
to the radial resolution is that it is hardly possible to de- 
scribe the protoneutron star atmosphere with about 100 
radial grid points in the context of a full supernova sim- 
ulation. In model CS we used the same 15 spatial grid 
points as in the Monte Carlo simulations. On the other 
hand, in model NI we used 105 spatial mesh points but 
the density, temperature and electron fraction were not 
interpolated between the grid points of the Monte Carlo 
simulations. 

In Fig. 16 we show the radial evolution of the aver- 
age energy as defined in Eq. (9) and that of the num- 
ber flux given by Eq. (19) for model CS, to be compared 
with the corresponding result for model ST2 in Fig. 14. 
It is clear that the agreement between the Monte Carlo 
and the Boltzmann results for both quantities is good. 
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Fig. 15. The energy spectrum of electron-type anti-neutrinos 
calculated with 18 energy mesh points (model FE, filled trian- 
gles) for background model Wl. The solid line shows the result 
of the Monte Carlo simulation, and the solid squares represent 
the result for the corresponding standard model ST2 with only 
12 energy grid points. 



Fig. 16. The radial evolution of the average energy of the 
flux (upper panel) and the number flux (lower panel) of the 
electron-type anti-neutrino for model CS where 15 radial mesh 
points are used and the background model is Wl. The filled 
triangles show the results of the Monte Carlo simulation, while 
the solid lines represent the results of the Boltzmann simula- 
tion. 



We note also that the angular distribution as well as the 
energy spectrum are hardly affected by this change of the 
spatial resolution. Model NI agrees with the Monte Carlo 
data nearly perfectly (except for the problems with the 
angular distribution discussed in Sect. 4.3) after averag- 
ing over spatial mesh zones in accordance to the way the 
Monte Carlo data represent the transport result. Since the 
Boltzmann results do not change with the number of radial 
grid points, we conclude that the quality of the numerical 
solutions is not degraded very much for simulations with 
a decreased spatial resolution. 

Minor oscillations of the number flux near the inner 
boundary can be seen in Fig. 14. This problem results 
from the fact that one cannot consistently specify the dis- 
tribution of neutrinos which leave the computational vol- 
ume at the inner boundary. While this distribution should 
be determined from the transport result just above the in- 
ner boundary the Boltzmann solver, however, requires an 
ad hoc specification in order to calculate the flux at the 
inner boundary. This leads inevitably to an inconsistency 
of the flux in the innermost zone and thus to the observed 
oscillations. In fact, when an inhomogeneous spatial mesh 
was used in model NU, in which the innermost grid zone 



was five times smaller than in the standard models, the 
oscillations were diminished as well. 

Finally, we illustrate possible errors which are associ- 
ated with the treatment of the finite differencing of the 
spatial advection term in the Boltzmann solver. In the 
radial advection term a linear average of the centered dif- 
ference and of the upwind difference is used with a weight 
factor that changes according to the ratio of the neu- 
trino mean free path to some chosen length scale. Mez- 
zacappa et al. (1993a,b,c) took the ratio of the mean free 
path to the local mesh width in order to construct the 
weighting. However, we found that this does not work well 
if the mesh width becomes of the same order as the mean 
free path but is much smaller than the scale height of the 
background. This is indeed the case in the inner optically 
thick region of our standard models with 105 radial zones. 
In a more recent version of his code, Mezzacappa (private 
communciation and 1998) defines the weighting factors by 
refering them to the neutrinosphcric radius. In Fig. 17 the 
dashed line shows the number flux of muon-type neutri- 
nos for model DI which used the prescription suggested by 
Mezzacappa et al. (1993a,b,c). The flux is slowly increas- 
ing with radius because the upwind differencing is given 
too large a contribution in the optically thick region where 



Shoichi Yamada et al.: Neutrino transport in type II supernovae: Boltzmann solver vs. Monte Carlo method 



17 



the centered differencing should actually be chosen. As 
demonstrated by the solid line in Fig. 17, the constancy 
of the flux, however, is recovered when the ratio of the 
mean free path to the distance up to the surface is chosen 
instead of the ratio of the mean free path to the local mesh 
width. Yet, this issue is probably not very important for 
realistic calculations of the entire neutron star, since the 
mesh width is usually not much smaller than the typical 
scale height of the matter distribution. 
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Fig. 17. The number fluxes of the muon-type neutrino 
for model DI (dashed line) and the corresponding standard 
model ST4 (solid line). The triangles are the result of the 
Monte Carlo simulation. 

To finish, we comment briefly on a last model in which 
we assumed that the nucleon scattering is taken isotropic 
to see to what extent the result changes. No significant 
effect was found by modifying the angular distribution of 
the dominant scattering reaction. 

5. Summary 

In this paper an extensive comparison was made between 
a newly developed Boltzmann neutrino transport code 
based on the discrete ordinate (Sat) method as described 
by Mezzacappa & Bruenn (1993a,b,c), and a Monte Carlo 
transport treatment (Janka 1987, 1991a) by performing 
time-dependent calculations of neutrino transport through 
realistic, static profiles of protoneutron star atmospheres 
under the assumption of spherical symmetry. In particu- 



lar, the sensitivity of the results of the Boltzmann solver 
to the employed numbers of radial, angular and energy 
grid points and to the treatment of the radial advection 
terms was investigated. The flux factor and Eddington fac- 
tor, which contain information about the angular distri- 
bution of the neutrinos in the neutrino-decoupling region, 
were also compared with the approximate treatment of 
this regime by a multi-energy-group flux-limited diffusion 
code (MGFLD). 

The Boltzmann and Monte Carlo results showed excel- 
lent agreement for observables such as the luminosity and 
the flux spectra which are determined in those regions of 
the star where the neutrino-matter interactions are still 
very frequent and thus the neutrino distributions are still 
nearly isotropic. Since the luminosity and the spectra are 
essentially conserved farther out, the spatial evolutions as 
well as the surface values exhibit this agreement, as long 
as the finite differencing of the Boltzmann solver is done 
in a conservative way. 

Some problems, however, were observed concerning the 
description of the angular distribution of the neutrinos by 
the Boltzmann results in the semitransparent and trans- 
parent regimes. Due to severe limitations of the number 
of angular grid points which can be used — typically only 
6-10 angular bins between zero and 180 degrees are com- 
patible with the steep increase of the requirements of com- 
puter time for better resolution — the Boltzmann code is 
not able to describe strong forward peaking of the neutrino 
distributions very well, if a quadrature set is employed for 
the angular integration where the maximum angle cosine 
Mmax is significantly less than unity. In this case the Boltz- 
mann results underestimate the anisotropy in the optically 
thin region outside the average "neutrinosphere" , and the 
exact limits for the flux factor, (fi) — > 1, and for the Ed- 
dington factor, (^i 2 ) — > 1, at large distances away from the 
neutrino source cannot be satisfactorily reproduced. This 
is in agreement with the trends also seen in recent results 
of Messer et al. (1998) and is exactly opposite to the defi- 
ciencies of MGFLD which tends to overestimate the radial 
beaming of the radiation because flux limiters enforce the 
free-streaming limit when the optical depth of the stellar 
background becomes very low (Janka 1991a, 1992). 

Since the dominant energy deposition rate by absorp- 
tion of electron neutrinos and antincutrinos in the hot- 
bubble region of the supernova core is inversely pro- 
portional to the flux factor (see Eq. (1)), which means 
that the energy transfer from neutrinos to the stellar 
plasma scales with the neutrino energy density (or num- 
ber density) in the heating region, one cannot exclude 
that the Boltzmann solver may lead to an overestima- 
tion of the neutrino heating in disadvantageous situa- 
tions, whereas MGFLD yields rates which are definitely 
too low. For a set of typical post-core bounce situations, 
Messer et al. (1998), however, claim on grounds of numer- 
ical tests that the net neutrino heating is converged with 
Sq and that the differences between and S4 are minor. 
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The problems may be more serious for neutrino reactions 
like neutrino-antineutrino annihilation which are sensitive 
to both the flux factor (/x) and the Eddington factor (/z 2 ) 
of the neutrinos (see Janka 1991b). 

The description of the angular neutrino distribution 
by the Boltzmann solver and thus the agreement with the 
highly accurate Monte Carlo data can be significantly im- 
proved by employing a variable angular mesh, even with- 
out increasing the total number of angular mesh points. 
The positions of the angular grid points must be moved 
at each time step and in every spatial zone such that they 
are clustering in the forward direction in the optically thin 
regime. 

We found that the energy spectra can be well calcu- 
lated with about 15 energy mesh points. The fact that a 
reduction of the number of spatial grid points from more 
than 100 to only 15 in the neutron star atmosphere did not 
change the quality of the transport results means that the 
Boltzmann code can be reliably applied to realistic simu- 
lations which involve the whole supernova core. Moreover, 
it was demonstrated that the details of the interpolation 
between centered differencing and upwind differencing in 
the spatial advection term can affect the accuracy of the 
transport results. 

The excellent overall agreement of the results obtained 
with the Boltzmann code and the Monte Carlo method 
confirms the reliability of both of them. Good performance 
of the S n method for solving the Boltzmann equation of 
neutrino transport has been found before by Mezzacappa 
& Bruenn (1993a,b,c) and Messer et al. (1998) even for 
realistic dynamic situations. We hope that the work de- 
scribed here also helps to reveal possible deficiencies and 
weaknesses and thus will serve for further improvements 
of the numerical treatment of neutrino transport in super- 
novae and protoneutron stars. 
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